*------------------------------------------------------------------------------
*
*By:  David E. Lewis
*Date:  July 5, 2000
*
*My hope is to read the file in to stata, take care of missing data
*set the data for survival analysis, and estimate all of the models for my
*myth of agency immortality paper.
*------------------------------------------------------------------------------
*First, I need to do some housekeeping.
set memory 64000
set l 72
set level 95
set matsize 150

*Now, I'll look at the spell file.
use c:\stata\durdata9.dta

*Now I'll describe the data.
describe

*Now I will take care of some missing cases.
keep if cumdur~=.
keep if year~=.
keep if estate~=.

*cumdur is cumulative duration or age.  year is what is appears to be.
*estate is end state.  It is 0,1 depending on whether an agency is active
*or terminated.

*summarize
 
*First, I need to take care of one problem which is that if an 
*agency was terminated on the first of the year its cumdur is 0 instead of
*1.  I also need to take care of the problem that if an agency was
*terminated on the first of the year that its cumdur=cumdur for the 
*previous case.
replace cumdur=1 if cumdur==0
replace cumdur=cumdur+1 if cumdur[_n-1]==cumdur

*Now I will try to get stata to recognize the spell file.  stset
*tell stata this is a spell file, cumdur is the agency age in days
*id is the unique agency id code.  There are multiple observations on single
*agencies.  Scale tells stata to report the results in years.

stset cumyear estate, id(agencyid) 

*Now I will describe the survival-time data.

stdes
more

*Now I will summarize the data

stsum
more

*I have to make the difference measure an absolute value.
generate abdistnc=abs(cshmed-cspres)
generate abhsecha=abs(hsechang)
generate abpresch=abs(preschan)

*Now with adjusted ADA scores.
generate adadist=abs(adahmed-presada)
generate abadhsc=abs(stadahm-adahmed)
generate abadprc=abs(stpres-presada)


*Now I have a problem with both measures of divergent preferences between the two branches that goof up the interactions.  That is, the larger
*the value, the lower the hazard rate.  To simplify things I take the maximum value of the distance measure and subtract the absolute value of the
*the difference in the scores.  This means an increase in the value (less preference divergence) should increase the hazard rate like the
*preference change measures.  From summarizing the data I know that the max of abdistnc=0.566 and adadist=37.

*So, in what follows I create the equivalent of the unified government measure for preferences.
generate nabdistnc=.566-abdistnc
generate nadadist=37-adadist


*Now I will change the party measures to republican indicators
generate reppres=1 if dempres==0
replace reppres=0 if dempres==1
generate rephouse=1 if demh==0
replace rephouse=0 if demh==1
*This next variable is the real growth in government agencies
generate realnum=counleg+counpma-terminat

*Finally an interaction for the party measures to show that the more dramatic the party turnover, the higher the hazard rate.
generate partintr=unfrmaj*unfrpres*unified


stcox unfrmaj unfrpres unified unemp war newadmin realnum temp reppres rephouse line leg, robust scaledsch(sca*) schoenfeld(sch*) nohr
stphtest, rank detail
predict outcome
summarize unemp realnum if outcome~=.

stcox unfrmaj unfrpres unified partintr unemp war newadmin realnum temp reppres rephouse line leg, robust nohr

more
stcox abhsecha abpresch abdistnc unemp war newadmin realnum temp cspres cshmed line leg, nohr
lrtest, saving(0)

*Now I will generate the interaction terms.  This is a little tricky since the amount of preference change should be in the same direction.

*Now for the preference measures.
generate hsecha=(bcshmed-cshmed)
generate prescha=(bcspres-cspres)

generate adhsc=(stadahm-adahmed)
generate adprc=(stpres-presada)

*Note: negative is getting more democratic and positive is more conservative.

generate prefinter=abhsecha*abpresch if hsecha>0&prescha>0|hsecha<0&prescha<0
replace prefinter=0 if hsecha>0&prescha<0|hsecha<0&prescha>0|hsecha==0|prescha==0

generate adainter=abadhsc*abadprc if adhsc>0&adprc>0|adhsc<0&adprc<0
replace adainter=0 if adhsc>0&adprc<0|adhsc<0&adprc>0|adhsc==0|adprc==0

stcox abhsecha abpresch nabdistnc prefinter unemp war newadmin realnum temp cspres cshmed line leg, nohr
more
lrtest, saving(1)
lrtest, using(0)
stcox abhsecha abpresch nabdistnc unemp war newadmin realnum temp cspres cshmed line leg, nohr
lrtest, using(1)
more
predict proutco 
summarize abhsecha abpresch nabdistnc unemp realnum cspres cshmed if proutco~=.

*Now with robust standard errors.
stcox abhsecha abpresch nabdistnc prefinter unemp war newadmin realnum temp cspres cshmed line leg, robust nohr

*Now with 1992 budget.
stcox unfrmaj unfrpres unified unemp war newadmin realnum temp reppres rephouse budget leg, robust  nohr
stcox unfrmaj unfrpres unified partintr unemp war newadmin realnum temp reppres rephouse budget leg, robust nohr 

*Now with adjusted ADA scores.
stcox nadadist abadhsc abadprc adainter unemp war newadmin realnum temp presada adahmed line leg, robust


*The following are the commands I used to estimate the graphs in the paper but the syntax seems sensitive to the version of stata being run.
*stcoxplt, xvar(unfrpres) adj(unfrmaj=0 partintr=0 unified=0 unemp war=0 newadmin=0 realnum temp=0 reppres=0 rephouse=0 line=1 leg=1) t1(Unfriendly President) l1(Cox Model Survival Probabilities) l2(" ") b2(Years) c(ll) s(ii) xlabel(0, 10, 20, 30, 40, 50)
*stcoxplt, xvar(unfrmaj) adj(partintr=0 unfrpres=0 unified=0 unemp war=0 newadmin=0 realnum temp=0 reppres=0 rephouse=0 line=1 leg=1) t1(Unfriendly Majority) l1(Cox Model Survival Probabilities) l2(" ") b2(Years) c(ll) s(ii) xlabel(0, 10, 20, 30, 40, 50)
*stcoxplt, xvar(unified) adj(unfrmaj=0 unfrpres=0 partintr=0 unemp war=0 newadmin=0 realnum temp=0 reppres=0 rephouse=0 line=1 leg=1) t1(Unified Government) l1(Cox Model Survival Probabilities) l2(" ") b2(Years) c(ll) s(ii) xlabel(0, 10, 20, 30, 40, 50)
*stcoxplt, xvar(partintr) adj(unfrmaj=1 unfrpres=1 unified=1 unemp war=0 newadmin=0 realnum temp=0 reppres=0 rephouse=0 line=1 leg=1) t1(Interaction of All 3 Factors) l1(Cox Model Survival Probabilities) l2(" ") b2(Years) c(ll) s(ii) xlabel(0, 10, 20, 30, 40, 50)






